function J = J_n( m, param)

    % Calculates J in natural wastage region

    J = (1-param.omega_1).*m./rho_fun( 1, 'N', param) - param.omega_0./rho_fun( 0, 'N', param) ...
            + param.J_1*m.^(param.gamma_N_1) + param.J_2*m.^(param.gamma_N_2);

end